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Abstract. We present a novel method for calculating the primordial non-Gaussianity 
produced by super-horizon evolution during inflation. Our method evolves the 
distribution of coarse-grained inflationary field values using a transport equation. 
We present simple evolution equations for the moments of this distribution, such as 
the variance and skewness. This method possesses some advantages over existing 
techniques. Among them, it cleanly separates multiple sources of primordial non- 
Gaussianity, and is computationally efficient when compared with popular alternatives, 
such as the 6N framework. We adduce numerical calculations demonstrating that our 
new method offers good agreement with those already in the literature. We focus on 
two fields and the /nl parameter, but we expect our method will generalize to multiple 
scalar fields and to moments of arbitrarily high order. We present our expressions in 
a field-space covariant form which we postulate to be valid for any number of fields. 

Keywords: Inflation, Cosmological perturbation theory. Physics of the early 
universe. Quantum field theory in curved spacetime. 
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1. Introduction 

Inflation generically predicts a primordial spectrum of density perturbations which 
is almost precisely Gaussian [H [2|, [3l HI [5]. In recent years the small non-Gaussian 
component [Gj [Tl [8], [HI HOl HH 112] has emerged as an important observable [13], and 
will be measured with good precision by the Planck Surveyor satellite [H] . In the near 
future, as observational data become more plentiful, it will be important to understand 
the non-Gaussian signal expected in a wide variety of models, and to anticipate what 
conclusions can be drawn about early-universe physics from a prospective detection of 
primordial non-Gaussianity. 

In this paper, we present a novel method for calculating the primordial non- 
Gaussianity produced by super-horizon evolution in two-field models of inflation. Our 
method is based on the real-space distribution of inflationary field values on a flat 
hypersurface, which can be thought of as a probability density function whose evolution 
is determined by a form of the coUisionless Boltzmann equation. Using a cumulant 
representation [HI [161 [T71 [H] to expand our density function around an exact Gaussian, 
we derive ordinary differential equations which evolve the moments of this distribution. 
Further, we show how these moments are related to observable quantities, such as the 
dimensionless bispectrum measured by /nl [H HH] ■ We present numerical results which 
show that this method gives good agreement with other techniques. It is not necessary to 
make any assumptions about the inflationary model beyond requiring a canonical kinetic 
term and applying the slow-roll approximation. While there are already numerous 
methods for computing the super- horizon contribution to /nl, including the widely 
used 5N formalism, we believe the one reported here has a number of advantages. 

First, this new technique is ideally suited to unraveling the various contributions 
to /nl- This is because we follow the moments of the inflaton distribution directly, 
which makes it straightforward to identify large contributions to the skewness or other 
moments. The evolution equation for each moment is simple and possesses clearly 
identifiable source terms, which are related to the properties of the inflationary flow 
on field space. This offers a clear separation between two key sources of primordial 
non-Gaussianity. One of these is the intrinsic non-linearity associated with evolution 
of the probability density function between successive flat hypersurfaces; the other is a 
gauge transformation from field fluctuations to the curvature peturbation, C,. It would 
be difficult or impossible to observe this split within the context of other calculational 
schemes, such as the conventional 5N formalism. 

A second advantage of our method is connected with the computational cost of 
numerical implementation. Analytic formulas for /nl are known in certain cases, mostly 
in the context of the 5N framework, but only for very specific choices of the potential 
[T9| [20l [2T| [22] or Hubble rate [23l [M] . These formulas become increasingly cumbersome 
as the number of fields increases, or if one studies higher moments [25| [26] . In the 
future, it seems clear that studies of complex models with many fields will increasingly 
rely on numerical methods. The numerical 5N framework requires the solution to a 
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number of ordinary differential equations which scales exponentially with the number 
of fields. Since some models include hundreds of fields, this may present a significant 
obstacle [27]. Moreover, the SN formalism depends crucially on a numerical integration 
algorithm with low noise properties, since finite differences must be extracted between 
perturbatively different initial conditions after ~ 60 e-folds of evolution. Thus, the 
background equations must be solved to great accuracy, since any small error has 
considerable scope to propagate. 

In this paper we ultimately solve our equations numerically to determine the 
evolution of moments in specific models. Our method requires the solution to a number 
of differential equations which scales at most polynomially (or in certain cases perhaps 
even linearly) with the number of fields. It does not rely on extracting finite differences, 
and therefore is much less susceptible to numerical noise. These advantages may be 
shared with other schemes, such as the numerical method recently employed by Lehners 
& Renaux-Petel [28] . 

A third advantage, to which we hope to return in a future publication, is that our 
formalism yields explicit evolution equations with source terms. From an analysis of 
these source terms, we hope that it will be possible to identify those physical features 
of specific models which lead to the production of large non-Gaussianities. 

This paper is organized as follows. In ^ we show how the non-Gaussian parameter 
/nl can be computed in our framework. The calculation remains in real space 
throughout (as opposed to Fourier space), which modifies the relationship between /nl 
and the multi-point functions of the inflaton field. Our expression for /nl shows a 
clean separation between different contributions to non-Gaussianity, especially between 
the intrinsic nonlinearity of the field evolution and the gauge transformation between 
comoving and flat hypersurfaces. In ^ we introduce our model for the distribution 
of inflaton field values, which is a "moment expansion" around a purely Gaussian 
distribution. We derive the equations which govern the evolution of the moments of 
this distribution in the one- and two-field cases. In §U we present a comparison of 
our new technique and those already in the literature. We compute /nl numerically in 
several two-field models, and find excellent agreement between techniques. We conclude 
in gi 

Throughout this paper, we use units in which c = h = 1, and the reduced Planck 
mass Mp ^ = SvrG is set to unity. 

2. Frameworks for computing /nl 

In this section, we introduce our new method for computing the non-Gaussianity 
parameter /nl- This method requires three main ingredients: a formula for the curvature 
perturbation, (, in terms of the field values on a spatially fiat hypersurface; expressions 
for the derivatives of the number of e-foldings, N, as a function of field values at horizon 
exit; and a prescription for evolving the field distribution from horizon exit to the time 
when we require the statistical properties of (. The first two ingredients are given in 
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Eqs. ©-(Uni) and found at the end of §221 and §231 respectively. The final 

ingredient is discussed in ^ 

2.1. Calculations beyond linear order 

Once it became clear that non-linearities of the microwave background anisotropics 
could be detected by the WMAP and Planck survey satellites [H], many authors 
studied higher-order correlations of the curvature perturbation. In early work, direct 
calculations of a correlation function were matched to the known limit of local non- 
gaussianity P, Uni EHl [301 EH [32]. This method works well if isocurvature modes are 
absent, so that the curvature perturbation is constant after horizon exit. In the more 
realistic situation that isocurvature modes cause evolution on superhorizon scales, all 
correlation functions become time dependent. Various formalisms have been employed 
to describe this evolution. Lyth & Rodriguez |Tl] extended the 6N method |33|, [M] 
beyond linear order. This method is simple and well-suited to analytical calculation. 
Rigopoulos, Shellard and van Tent [351 [3S] worked with a gradient expansion, rewriting 
the field equations in Langevin form. The noise term was used as a proxy for setting 
initial conditions at horizon crossing. A similar 'exact' gradient formalism was written 
down by Langlois & Vernizzi [37l|3Hl[39j. In its perturbative form, this approach has been 
used by Lehners & Renaux-Petel to obtain numerical results [28]. Another numerical 
scheme has been introduced by Huston & Malik |40j . 

What properties do we require of a successful prediction? Consider a typical 
observer, drawn at random from an ensemble of realizations of inflation. In any of 
the formalisms discussed above, we aim to estimate the statistical properties of the 
curvature perturbation which would be measured by such an observer. Some realizations 
may yield statistical properties which are quite different from the ensemble average, but 
these large excursions are uninteresting unless anthropic arguments are in play. 

Next we introduce a collection of comparably-sized spacetime volumes whose 
mutual scatter is destined to dominate the microwave background anisotropy on a given 
scale. Neglecting spatial gradients, each spacetime volume will follow a trajectory in field 
space which is slightly displaced from its neighbors. The scatter between trajectories is 
determined by initial conditions set at horizon exit, which are determined by promoting 
the vacuum fluctuation to a classical pcrturbationg A correct prediction is a function 
of the trajectories followed by every volume in the collection, taken as a whole. One 
never makes a prediction for a single trajectory. 

Each spacetime volume follows a trajectory, which we label with its position 

f The scatter among trajectories may grow with the overall volume, owing to back-reaction effects in 
de Sitter space [4T1 l42l l43l |44] . If inflation can end in more than one vacuum, leading to different 
predictions for observable quantities, it may be helpful to evaluate this enhanced dispersion for the 
purpose of deciding into which vacuum the field will fall. Once this minimum has been decided, however, 
large-volume effects needlessly complicate the calculation. In this paper, we will always assume that 
predictions are being made by following a sufficient number of trajectories to determine the statistical 
properties with reasonable precision, but no more. 
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at some fixed time, to be made precise below. Throughout this paper, superscript 
denotes evaluation on a spatially flat hypersurface. Consider the evolution of some 
quantity of interest, F, which is a function of trajectory. If we know the distribution 
P{ip*) we can study statistical properties of F such as the m^^ moment Km, 



p{^^) [F{^*) - {F)r , (1) 

where we have introduced the ensemble average of F, 

/ d^^P(^^)F(^^). (2) 



In Eqs. ([I])-(l2]), if* stands for a collection of any number of fields. It is the Km which 
are observable quantities. 

Eq. ([T]) defines what we will call the exact separate universe picture. It is often 
convenient to expand F{ip*) as a power series in the field values centered on a fiducial 
trajectory, labelled 'fid,' 



(3) 



n=l 

When Eq. ([3]) is used to evaluate the Km, we refer to the 'perturbative' separate universe 
picture. If all terms in the power series are retained, these two versions of the calculation 
are formally equivalent. In unfavorable cases, however, convergence may occur slowly 
or not at all. This possibility was discussed in Refs. [121 US]- Although our calculation 
is formally perturbative, it is not directly equivalent to Eq. ([3]). We briefly discuss the 
relation of our calculation to conventional perturbation theory in ^ 



2.2. Calculating /nl in the separate universe picture 

By definition, the curvature perturbation ( measures local fluctuations in expansion 
history (expressed in e- folds N), calculated on a comoving hypersurface. In many 
models, the curvature perturbation is synthesized by superhorizon physics, which 
reprocesses a set of Gaussian fluctuations generated at horizon exit. In a single-field 
model, only one Gaussian fluctuation can be present, which we label (g. Neglecting 
spatial gradients, the total curvature perturbation must then be a function of (g alone. 
For \(g\ <^ 1, this can be well-approximated by 

C =^ C. + ^/nl (C,^ - (C,^)) , (4) 

where /nl is independent of spatial position. Eq. (jl]) defines the so-called "local" form 
of non-gaussianity. It applies only when quantum interference effects can be neglected, 
making ( a well-defined object rather than a superposition of operators ^7\. If this 
condition is satisfied, spatial correlations of ( may be extracted and it follows that /nl 
can be estimated using the rule 

^ 18 W 
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where we have recalled that ( is nearly Gaussian, or equivalently that |/nl| ^ lCgl~^- 

With /nl spatially independent, Eq. (jl]) strictly applies only in single-field infiation. 
In this case one can accurately determine /nl by applying Eq. (jl]) to a single trajectory 
with fixed initial conditions, as in the method of Lehners & Renaux-Petel |28j. Where 
more than one field is present, /nl may vary in space because it depends on the 
isocurvature modes. In this case one must determine /nl statistically on a bundle of 
adjacent trajectories which sample the local distribution of isocurvature modes. Eq. ([5]) 
is then indispensible. Following Maldacena ^Q\, and later Lyth & Rodriguez [11], we 
adopt Eq. ([5]) as our definition of /nl, whatever its origin. In real space, the coefficient 
5/18 in Eq. ([5]) depends on the convention {Q = 0. More generally, this follows from the 
definition of k^, Eq. ([1]). In Fourier space, either prescription is automatically enforced 
after dropping disconnected contributions, again leading to Eq. ([5]). 

To proceed, we require estimates of the correlation functions {(() and (CCO- We 
first describe the conventional approach, in which '^ir' denotes a fiat hypersurface at a 
fixed initial time. The quantity N{ip*) denotes the number of e-foldings between this 
initial slice and a final comoving hypersurface, where i indexes the species of light scalar 
fields. The local variation in expansion can be written in the fiducial picture as 



oo ^ 



6vl{^)---6vl{^), (6) 

where Sif* = if* - (/^^gj. 

Subject to the condition that the relevant scales are all outside the horizon, we are 
free to choose the initial time — set by the hypersurface — at our convenience. In the 
conventional approach, is taken to lie a few e-folds after our collection of spacetime 
volumes passes outside the causal horizon [TTIIT2] . This choice has many virtues. First, 
we need to know statistical properties of the field fiuctuations 6ip* only around the 
time of horizon crossing, where they can be computed without the appearance of large 
logarithms |^H8]. Second, as a consequence of the slow-roll approximation, the 6ip* are 
uncorrelated at this time, leading to algebraic simplifications. Finally, the 6N formula 
subsumes a gauge transformation from the field variables 6(f* to the observational 
variable (. Using Eqs. ([I])-(l2l), ([5]) and ([6]), one finds that /nl can be written to a 
good approximation [11] 



/nl ~ ^ ... xo ) * at horizon crossing (7) 

6 {N^kN^kV 



where iV j = dN/ dip* and for simplicity we have dropped the ^-k^ which indicates time of 
evaluation. A similar definition applies for N^ij. 

Eq. (I7j) is accurate up to small intrinsic non-Gaussianities present in the field 
fiuctuations at horizon exit. As a means of predicting /nl it is pleasingly compact, 
and straightforward to evaluate in many models. Unfortunately, it also obscures the 
physics which determines |/nl|- For this reason it is hard to infer, from Eq. ([7]) alone. 
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those classes of models in which |/nl| is always large or small||| 

Our strategy is quite different. We choose to lie around the time when we require 
the statistical properties of C,. The role of the 5N formula, Eq. ([6]), is then to encode 
only the gauge transformation between the 5^* and C,. In §2.31 below, we show how the 
appropriate gauge transformation is computed using the 5N formula. In the present 
section we restrict our attention to determining a formula for /nl under the assumption 
that the distribution of field values on '^^r' is known. In [|3l we will supply the required 
prescription to evolve the distribution of field values between horizon exit and ^-k\ 

Although the 5ip\ are independent random variables at horizon exit, correlations 
can be induced by subsequent evolution. One must therefore allow for off-diagonal 
terms in the two-point function. Remembering that we are working with a collection of 
spacetime volumes in real space, smoothed on some characteristic scale, we write 

{5^*,5v)) ^ S,,. (8) 

Sjj does not vary in space, but it may be a function of the scale which characterizes our 
ensemble of spacetime volumes. In all but the simplest models it will vary in time. It 
is also necessary to account for intrinsic non-linearities among the which are small 
at horizon crossing but may grow. We define 

{5^15^)5^1) = a,jk. (9) 

Likewise, aijk should be regarded as a function of time and scale. The permutation 
symmetries of an expectation value such as ([9]) guarantee that, for example, ai22 = 
^212 = a22ili When written explicitly, we place the indices of symbols such as a in 
numerical order. Neglecting a small (< 0(S^)) intrinsic four-point correlation, it follows 
that 

Now we specialize to a two-field model, parametrized by fields (pi and ip2- Using 
Eqs. (II])-(l2l), (ini) and ([H]), it follows that the two-point function of C satisfies 

(CO = A^'Sii + NlT.22 + 2NiN2^i2 arbitrary (11) 

The three-point function can be written 

(CCC) = (CCC)i + (CCC)2, (12) 

where we have identified two separate contributions, labelled '1' and '2'. The '1' term 
includes all contributions involving intrinsic non-linearities, those which arise from non- 
Gaussian correlations among the field fiuctuations, 

(CCC)i = + iV,|a222 + 3NlN 2au2 + 3NiNlau2. arbitrary (13) 

I Even in simple models it can be quite subtle to determine what range of |/nl| is dynamically allowed. 
See, for example, Refs. [2T|[22]. 

§ We are assuming that these expectation values are ensemble averages over classical stochastic fields, 
and are therefore invariant under reordering of the fields. 
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The '2' term encodes non-linearities arising directly from the gauge transformation to C, 
(CCC)2 = ^iV_iiV,nS2, + ^-NlN,22T.l^ W arbitrary 

+ 9 {N^iN^2N,ii + NlN^i2) S11S12 + 9 (AT, 1 AT, 2 A^, 22 + A^^ Ar,i2) S12S22 

+ ^ (AT^iV^n + Ar2iV,22 + 4Ar,iAr,2Ar,i2) (SnS22 + 

(AT^nSn + 2Ar,i2Si2 + Ar,22S22) (CO, (14) 

After use of Eq. ([5]), Eqs. (|T2l) - (|T^ can be used to extract the non-linearity parameter 
/nl- This decomposes likewise into two contributions /nl = /nli + /nl2, which we shall 
discuss in more detail in §11 

2.3. The derivatives of N 

To compute /nl in concrete models, we require expressions for the derivatives N^i and 
N^ij. For generic initial and final times, these are difficult to obtain. Lyth & Rodriguez 
[TT] used direct integration, which is effective for quadratic potentials and constant slow- 
roll parameters. Vernizzi & Wands [I9] obtained expressions in a two-field model with 
an arbitrary sum-separable potential by introducing Gaussian normal coordinates on 
the space of trajectories. Their approach was generalized to many fields by Battefeld 
& Easther [20]. Product-separable potentials can be accommodated using the same 
technique |49]. An alternative technique has been proposed by Yokoyama et al. [50]. 

A considerable simplification occurs in the present case, because we only require 
the derivative evaluated between fiat and comoving hypersurfaces which coincide in 
the unperturbed universe. For any species i, and to leading order in the slow-roll 
approximation, the number of e-folds A^ between the fiat hypersurface and a comoving 
hypersurface 'c' satisfies 

N = — — dipi no sum on i, (15) 

where Vi = dV/dipi and {{p*,{p'^} are the field values evaluated on and 'c,' 
respectively. Under an infinitesimal shift of (p*, we deduce that A^^j obeys 

"■'-{vj-ivji^ <^^) 

Note that this applies for an arbitrary V, which need not factorize into a sum or product 
of potentials for the individual species i. In principle a contribution from variation of 
the integrand is present, which spoils a naive attempt to generalize the method of 
Refs. [T9l [20t H9] to an arbitrary potential. This contribution vanishes in virtue of our 
supposition that and 'c' are infinitesimally separated. 

To compute dip'^/dip* it is helpful to introduce a quantity C, which in the sum- 
separable case coincides with the conserved quantity of Vernizzi & Wands [191 EI]. For 
our specific choice of a two-field model, this takes the form 



Ci^i,^2)^ I y;M- j y-^M, (17) 
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where the integrals are evaluated on a single spatial hypersurface. In an M-field 
model, one would obtain M — 1 conserved quantities which label the isocurvature fields. 
The construction of these quantities is discussed in Refs. [Sni [25] . For sum-separable 
potentials one can show using the equations of motion that C is conserved under time 
evolution to leading order in slow-roll. It is not conserved for general potentials, but 
the variation can be neglected for infinitesimally separated hyper surf aces. 
Under a change of trajectory, C varies according to the rules 

l^ = - (18) 

and 

dC _ _H^ 
dip2 V2 ' 
The comoving hypersurface 'c' is defined by 

- + 1^1) +V = constant. 



(19) 



(20) 



We are assuming that the slow-roll approximation applies, so that the kinetic energy 
may be neglected in comparison with the potential V. Therefore on 

0. 



dV difl ^ dV d^l 



di^l dC di^l dC 

Combining Eqs. (fT8|) . (fT9l) and (!2T!) we obtain expressions for difl/dip'*, namely 



dip* 

where we have defined 
tan^ e = 



V 



v_ 



sin^ 



V 
V 



V 
V 



sin' 9, 



^\ {^ \ cos'9 



V 



iV,2 



V 

V^2 



COS 6, 



(21) 

(22) 
(23) 
(24) 
(25) 

(26) 



Eqs. fl22|) - fl25|) can alternatively be derived without use of C by comparing Eq. f|T6|) with 
the formulas of Ref. [52], which were derived using conventional perturbation theory. 
Applying (fT6l) . we obtain 



V 



COS 9\ 



V 

V^2 



To proceed, we require the second derivatives of A^. 
from ([27D, after use of Eqs. We find 

*2 



sin=^ e. (27) 
These can be obtained directly 



N 



11 



1 - 



X 



VVm 



cos' 9 + 2 



V 



sin^9 



V 

V,lV,2 



cos' 9 



Wo 



sin* 9 - 



V 



11 



22 



V 



V 



+ 



V,2V,12 

VVi 



COS 6* sin 9 



{21 
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An analogous expression for N^22 can be obtained after the simultaneous exchange 
{1^2, sin ^ cos}. The mixed derivative satisfies 



N 



12 



v_ 



v_ 



cos^^ 



sin^^ 



12 



+ cos^ Q 



V 



VV2 



sin^ e + 



V 



11 



V 



K22 _l_ K2K12 



V 



VVi 



COS 6' sin 9 



(29) 



Now that the calculation is complete, we can drop the superscripts and 'c,' since any 
background quantity is the same on either hypersurface. Once this is done it can be 
verified that (despite appearances) Eq. (1291) is invariant under the exchange 1 <r->- 2. 



3. Transport equations 

In this section we return to the problem of evolution between horizon exit and the 
time of observation, and supply the prescription which connects the distribution of field 
values at these two times. 



3.1. Non-gaussian distribution in the single-field case 

We begin by discussing the single-field system, which lacks the technical complexity of 
the two-field case, yet still exhibits certain interesting features which recur there. Among 
these features are the subtle difference between motion of the statistical mean and the 
background field value, and the hierarchy of moment evolution equations. Moreover, 
the structure of the moment mixing equations is similar to that which obtains in the 
two-field case. For this reason, the one-field scenario provides an instructive example of 
the techniques we wish to employ. 

Recall that we work in real space with a collection of comparably sized spacetime 
volumes, each with a slightly different expansion history, and the scatter in these 
histories determines the microwave background anisotropy on a given angular scale. 
Within each volume the smoothed background field ip takes a uniform value described 
by a density function P{'^), where in this section we are dropping the superscript '-k^ 
denoting evaluation of spatially fiat hypersurfaces. Our ultimate goal is to calculate the 
reduced bispectrum, /nl, which describes the third moment of P{^p). In the language 
of probability this is the skewness, which we denote a. A Gaussian distribution has 
skewness zero, and inflation usually predicts that the skew is small. For this reason, 
rather than seek a distribution with non-zero third moment, as proposed in Ref. |18j . 
we will introduce higher moments as perturbative corrections to the Gaussian. Such a 
procedure is known as a cumulant expansion. 

The construction of cumulant expansions is a classical problem in probability theory. 
We seek a distribution with centroid ip^, variance o"^, and skew a, with all higher 
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moments determined by cr and a alone. A distribution with suitable properties is 



cr 



where 



2a2 



(30) 



(31) 



is a pure Gaussian and Hn denotes the n*^ Hermite polynomial, for which there are 
multiple normalization conventions. We choose to normalize so that 

/oo -| 
-=e-^'/2 H^^x)HUx) dx = n\ (32) 

which implies that the leading term of Hn{x) is x". This is sometimes called the 
"Probabilist's convention." We define expectation values (■ ■ •) by the usual rule, 

/oo 
P((^)F dx. (33) 
-oo 

The probabihty density function in Eq. (l30ll has the propertied 

(1) = 1, {ip) = ipo, {{if - ifof) = a\ and {{if - ip,f) = a. (34) 

The moments y^o, a, and a may be time- dependent, so evolution of the probability 
density in time can be accommodated by finding evolution equations for these quantities. 

The density function given in Eq. ( l30l) is well-known and has been applied in 
many situations. It is a solution to the problem of approximating a nearly-Gaussian 
distribution whose moments are known. Eq. (l30l) is in fact the first two terms of 
the Gram-Charlier 'A ' series, also sometimes called the Gauss-Hermite series ^ In 
recent years it has found multiple applications to cosmology, of which our method 
is closest to that of Taylor & Watts [53]. Other applications are discussed in 
Refs. [ISlIITlIISllMlIlHlESlESlEHlETlES For a review of the 'A' series and 

related nearly-Gaussian probability distributions from an astrophysical perspective, see 
|61j . In this paper, we will refer to Eq. (I30l) and its natural generalization to higher 
moments as the "moment expansion." 

In the slow-roll approximation, the field in each spacetime volume obeys a simple 
equation of motion 

d^ d\nV{^) 



dN dip 



u{ip), (35) 



f These formulas apply for arbitrary values of a, and do not depend on the approximation that a is 
small. However, for large a the density function pop may become negative for some values of tp. It then 
ceases to be a probability density in the strict sense. This does not present a problem in practice, since 
we are interested in distributions which are approximately Gaussian, and for which a will typically be 
small. Moreover, our principal use of Eq. (|30|) is as a formal tool to extract evolution equations for 
each moment. For this reason we will not worry whether P{ip) defines an honest probability density 
function in the strict mathematical sense. 

I In the physics literature, this series has sometimes erroneously been called the Edgeworth expansion. 
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where N records the number of e- foldings of expansion. We refer to u{!f) as the velocity 
field. Expanding u about the instantaneous centroid ipo gives 

u{ip) =uo + u^iip - (po) + ^Mw(v5 - <Pof H , (36) 



where 



du 

dip 



d^M 



(37) 



The value of ipo evolves with time, so each expansion coefficient is time- dependent. 
Hence, we do not assume that the velocity field is globally well- described by a quadratic 
Taylor expansion, but merely that it is well-described as such in the neighborhood of 
the instantaneous centroid. We expand the velocity field to second order, although in 
principle this expansion could be carried to arbitrary order. 

It remains to specify how the probability density evolves in time. Conservation of 
probability leads to the transport equation 
dP d 

Eq. fl38|) can also be understood as the limit of a Chapman-Kolmogorov process as 
the size of each hop goes to zero. It is well known — for example, from the study 
of Starobinsky's diffusion equation which forms the basis of the stochastic approach 
to infiation [621 — that the choice of time variable in this equation is significant, with 
different choices corresponding to the selection of a temporal gauge. We have chosen 
to use the e-folding time, A^, which means that we are evolving the distribution on 
hypersurfaces of uniform expansion. These are the spatially fiat hypersurfaces whose 
field perturbations enter the SN formulas described in ^ 

In principle, Eq. (1551) can be solved directly. In practice it is simpler to extract 
equations for the moments of P, giving evolution equations for ipQ, a and a. To achieve 
this, one need only resolve Eq. fl38|) into a Hermite series of the form 

PgJ2cnHn{^-y^o) = (39) 

The Hermite polynomials are linearly independent, and application of the orthogonality 
condition (152]) shows that the c„ must all vanish. This leads to a hierarchy of equations 
Cn = 0, which we refer to as the moment hierarchy. At the top of the hierarchy, the 
equation cq = is empty and expresses conservation of probability. 

The first non-trivial equation requires Ci = and yields an evolution equation for 
the centroid ipo, 



dipo 1 2 



uq + -Mw*^^- (40) 



The first term on the right-hand side drives the centroid along the velocity field, as one 
would anticipate based on the background equation of motion, Eq. (I5S]) . However, the 
second term shows that the centroid is also infiuenced as the wings of the probability 
distribution probe the nearby velocity field. This infiuence is not captured by the 
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background equation of motion. If we are in a situation with u^p^p > 0, then the wings 
of the density function will be moving faster than the center. Hence, the velocity of the 
centroid will be larger than one might expect by restricting attention to (po- Accordingly, 
the mean fluctuation value is not following a solution to the background equations of 
motion. 

Evolution equations for the variance and skew a are obtained after enforcing 
C2 = C3 = 0, yielding 



2u^a +u^^a (41) 



3u^a + Su^^cr^ (42) 



dA^ 
da 

diV 

In both equations, the first term on the right-hand sides describes how a and a scale as 
the density function expands or contracts in response to the velocity field. These terms 
force 0"^ and a to scale in proportion to the velocity field. Specifically, if we temporarily 
drop the second terms in each equation above, one finds that ~ and a ~ u^. This 
precisely matches our expectation for the scaling of these quantities. Hence, these terms 
account for the Jacobians associated with infinitesimal transformations induced by the 
flow u{{p). 

For applications to inflationary non-Gaussianity, the second terms in fj4Tl) and P2l) 
are more relevant. These terms describe how each moment is sourced by higher moments 
and the interaction of the density function with the velocity field. In the example above, 
if we are in a situation where > 0, the tails of the density function are moving faster 
than the core. This means that one tail is shrinking and the other is extending, skewing 
the probability density. The opposite occurs when < 0. These effects are measured 
by the second term in ( 142|) . Hence, by expanding our PDF to the third moment, and our 
velocity field to quadratic order, we are able to construct a set of evolution equations 
which include the leading-order source terms for each moment. 

3.2. The two-field case 

There is little conceptually new as we move from one field to two. The new features are 
mostly technical in nature. Our primary challenge is a generalization of the moment 
expansion to two fields, allowing for the possibility of correlation between the fields. 
With this done, we can write down evolution equations whose structure is very similar 
to those found in the single-field case. 

The two-field system is described by a two-dimensional velocity field Ui, defined by 

= ^, (43) 

where again we are using the number of e-folds as the time variable. The index i 
takes values in {1,2}. While we think it is likely that our equations generalize to any 
number of fields, we have only explicitly constructed them for a two-field system. As 
will become clear below, certain steps in this construction apply only for two fields, and 
hence we make no claims at present concerning examples with three or more fields. 
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The two-dimensional transport equation is 

^,_|_,^.P(,..)H0. (44, 

Here and in the following we have returned to our convention that repeated species 
indices are summed. As in the single-field case, we construct a probability distribution 
which is nearly Gaussian, but has a small non-zero skewness. That gives 

P{ipi, N) = Pgiipi, N)Png{ipi, N) (45) 

where Pg is a pure Gaussian distribution, defined by 



27r Vdet E 



exp 



(46) 



In this equation, $j defines the center of the distribution and S describes the covariance 
between the fields. We adopt a conventional parametrization in terms of variances af 
and a correlation coefficient p, 



pa 1(72 

The matrix a defines two-point correlations of the fields. 



(47) 



(4^ 



All skewnesses are encoded in Png- Before defining this explicitly, it is helpful 



to pause and notice a complication inherent in Eqs. fHB]) - fH7|) which was not present 
in the single-field case. To extract a hierarchy of moment evolution equations from 
the transport equation, Eq. (135]) . we made the expansion given in (15^ and argued 
that orthogonality of the Hermite polynomials implied the hierarchy c„ = 0. However, 
Hermite polynomials of the form if„[((y9j — $j)/o"] are not orthogonal under the Gaussian 
measure of Eq. (H6ll . Following an expansion analogous to Eq. (!39l) the moment hierarchy 
would comprise linear combinations of the coefficients. The problem is essentially an 
algebraic question of Gram-Schmidt orthogonalization. 

To avoid this problem it is convenient to diagonalize the covariance matrix S, 
introducing new variables X and Y for which Eq. ( H6l) factorizes into the product of two 
measures under which the polynomials Hn{X) and Hn{Y) are separately orthogonal. 
The necessary redefinitions are 



and 



X 



Y = 



1 



V2(l + p) 



0"! 



(p2 - $2 



0-2 



%1-p) LV ^1 



^2 



0-2 



A simple expression for Pg can be given in terms of X and Y, 



(49) 



(50) 
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We now define the non-Gaussian factor, wliicli encodes tlie skewnesses, to be 

+ ^H,{X)H,{Y) + ^H.iXmY) + ^H,iY). (52) 
In tliese variables we find (X^) = (V^) = 1, but (XY) = 0. In addition, we have 
{XXX) = axxx, {XXY) = axxY, (XYY) = axYY, and (FFF) = «yyy. (53) 

In order for Eq. fl52l) to be useful, it is necessary to express the skewnesses associated 
with the physical variables ipi in terms of X and Y. By definition, these satisfy 

{{ifii - - $j)(<^fc - $fc)) = aijk. (54) 

After substituting for the definition of these quantities inside the expectation values in 
Eq. fl53p we arrive at the relations 

1 / Qui o ail2 o «122 



1 / am aii2 0122 "222 \ /^^x 

2v/2(l-p)(l + p) V 0^1 ^1(^2 ol ) 

1 / Q^lll _ Q112 _ 0:122 0222\ 



axyy = — , , — o n 2 ^ ' (5"^) 

2v/2(l +p)(l - p) V ^1 ^?^2 ^1^2 c^i / 

1 /am o '^112 , o'^122 0222 A /-ON 

oyyy = — 1= 7- — 5 i—^ V i ^ ^ • (58) 

2^2(1 - p)3/2 V (T? a\a2 a^ol o\ ) ^ ' 

The moments $j, Sjj and ajj^ are time-dependent, but for clarity we will usually 
suppress this in our notation. 

Next we must extract the moment hierarchy, which governs evolution of $i, cXi, p 
and ajjfc- We expand the velocity field in a neighborhood of the instantaneous centroid 
$j according to 

Ui{^j) = Uio + Uiji^Lpj - $j) + ^Uijki^Pj - ^j){^Pk - $fc) H , (59) 

where we have defined 

dui 



dip. 



and Uijk 



(60) 



As in the single-field case, these coefficients are functions of time and vary with the 
motion of the centroid. The expansion can be pursued to higher order if desired. 

Our construction of X and Y implies that the two-field transport equation can be 
arranged double Gauss-Hermite expansion, 

^^^^^ + ^ [Ui P(^„ N)] =PgY. CmnHm{X)H4Y) = 0. (61) 

m,n>0 

Because the Hermite polynomials are orthogonal in the measure defined by Pg, we 
deduce the moment hierarchy 

Cmn = 0. (62) 
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We define the "rank" r of eacli coefficient Cmn by r = m + n. We terminated the 
velocity field expansion at quadratic order, and our probability distribution included 
only the first three moments. It follows that only Cmn with rank five or less are nonzero. 
If we followed the velocity field to higher order, or included higher terms in the moment 
expansion, we would obtain non-trivial higher-rank coefficients. Inclusion of additional 
coefficients requires no qualitative modification of our analysis and can be incorporated 
in the scheme we describe below. 

A useful feature of the expansion in Eq. f l6Tl) is that the rank-r coefficients give 
evolution equations for the order-r moments. Written explicitly in components, the 
expressions that result from (16T1) are quite cumbersome. However, when written as 
field-space covariant expressions they can be expressed in a surprisingly compact form. 

Rank The rank-0 coefficient Cqo is identically zero. This expresses the fact that the 
total probability is conserved as the distribution evolves. 

Rank 1 The rank-1 coefficients Cqi and Cio give evolution equations for the centroid $j. 
These equations can be written in the form 

= Uio + -l^jkUijk- (63) 

We remind the reader that here and below, terms like Uio, Uij and Uijk represent 
the velocity field and its derivatives evaluated at the centroid <I>i. The first term in 
(l63l) expresses the non-anomalous motion of the centroid, which coincides with the 
background velocity field of Eq. (H3|) . The second term describes how the wings 
of the probability distribution sample the velocity field at nearby points. Narrow 
probability distributions have small components of S and hence are only sensitive 
to the local value of Ui{(fj). Broad probability distributions have large components 
of S and are therefore more sensitive to the velocity field far from the centroid. 

Rank 2 The rank-2 coefficients Co2, Cn and C20 give evolution equations for the variances 
af and the correlation p. These can conveniently be packaged as evolution equations 
for the matrix S 

dX] 1 

This equation describes the stretching and rotation of S as it is transported by the 
velocity field. It includes a sensitivity to the wings of the probability distribution, 
in a manner analogous to the similar term appearing in ( 1411) . Hence the skew aijk 
acts as a source for the correlation matrix. 

Rank 3 The rank-3 coefficients cqs, C12, C21 and C30 describe evolution of the moments 
aijk- These are 

+ cyclic permutations i j k. (65) 

The first term describes how the moments flow into each other as the velocity field 
rotates and shears the (X, Y) coordinate frame relative to the <fi coordinate frame. 
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The second term describes sourcing of non-Gaussianity from inhomogeneities in the 
velocity field and the overall spread of the probability distribution. 

Some higher-rank coefficients — in our case, those of ranks four and five — are also 
nonzero, but do not give any new evolution equations. These coefficients measure the 
"error" introduced by truncating the moment expansion. If we had included higher 
cumulants, these higher-rank coefficients would have given evolution equations for the 
higher moments of the probability distribution. In general, all moments of the density 
function will mix so it is always necessary to terminate our expansion at a predetermined 
order — both in cumulants and powers of the field fluctuation. The order we have chosen 
is sufficient to generate evolution equations containing both the leading-order behavior 
of the moments — namely, the first terms in Eqs. (!63l) . (!64l) and (!65l) — and the leading 
corrections, given by the latter terms in these equations. 

4. Numerical results 

At this point we put our new method into practice. We study two models for which the 
non-Gaussian signal is already known, using the standard 6N formula. For each case 
we employ our method and compare it with results obtained using 6N. To ensure a fair 
comparison, we solve numerically in both cases. Our new method employs the slow-roll 
approximation, as described above. Therefore, when using the SN approach we produce 
results both with and without slow-roll simplifications. 

First consider double quadratic inflation, which was studied by Rigopoulos, Shellard 
& van Tent [36l |63] and later by Vernizzi & Wands [19] . The potential is 

We use the initial conditions chosen in Ref. where m^/m^ = 9, and the flducial 
trajectory has coordinates = 8.2 and x* = 12.9. 

We plot the evolution of /nl in Fig. [1], which also shows the prediction of the 
standard 6N formula (with and without employing slow roll simpliflcations) . We 
implement the 6N algorithm using a flnite difference method to calculate the derivatives 
of A^. A similar technique was used in Ref. [19]. This model yields a very modest non- 
Gaussian signal, below unity even at its peak. If inflation ends away from the spike then 
/nl is practically negligible. 

Eq. (fT2l) shows that the method of moment transport allows us to separate 
contributions to /nl from the intrinsic non-Gaussianity of the fleld fluctuations, and 
non-linearities of the gauge transformation to As explained in §2.2[ we denote the 
former /nli and the latter /nl2, and plot them separately in Fig. [2l Inspection of 
this flgure clearly shows that /nl is determined by a cancellation between two much 
larger components. Its flnal shape and magnitude are exquisitely sensitive to their 
relative phase. Initially, the magnitudes of /nli and /nl2 grow, but their sum remains 
small. The peak in Fig. [1] arises from the peak of /nl2; which is incompletely cancelled 
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0.25 




Figure 1. Evolution of /nl in double quadratic inflation. The solid red line is obtained 
by numerically solving the moment transport equations obtained in Sj3l The blue 
dashed line and green dot-dashed line are the output of a numerical implementation of 
the standard 5N approach, with and without slow roll respectively, using the fiducial 
picture. The red and blue lines lie on top of each other. 
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Figure 2. Evolution of /nli (solid red line), measuring the contribution of intrinsic 
non-linearities among the field fluctuations; and /nl2 (dashed blue line), measuring 
the contribution of the gauge transformation to Q. 
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Figure 3. Evolution of /nl, for tiie potential V — Vox^e"^* (Example A from §5 
of Ref. i21]). The solid red line represents the method of moment transport, whereas 
the blue dashed line and green dot-dashed line represents the output of conventional 
numerical 5N with and without slow-roll respectively. The red and blue lines are again 
coincident. 



by /nli- It is remarkable that /nli initially evolves in exact opposition to the gauge 
transformation, to which it is not obviously connected. 

In the double quadratic model, /nl is always small. However, it has recently 
been shown by Byrnes et al. that a large non-Gaussian signal can be generated even 
when slow- roll is a good approximation [211 [22] . The conditions for this to occur are 
incompletely understood, but apparently require a specific choice of potential and strong 
tuning of initial conditions. In Figs. [3H1] we show the evolution of /nl in a model with 
the potential 

V = V,x^e-^'^\ (67) 

which corresponds to Example A of Ref. [21] §5] when we choose A = 0.05 and initial 
conditions x* = 16, 0* = 0.001. It is clear that the agreement is exact. In this 
model, /nl is overwhelmingly dominated by the contribution from the second-order 
gauge transformation, /nl2, as shown in Fig. HI This conclusion applies equally to the 
other large- /nl examples discussed in Refs. [2T], [22], although we make no claim that 
this is a general phenomenon. 

In conclusion. Figs. [1] and [3] show excellent agreement between our new method 
and the outcome of the numerical 5N formula. These figures also compare the moment 
transport method and 6N without the slow-roll approximation. We conclude that the 
slow-roll estimate remains broadly accurate throughout the entire evolution. 
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Figure 4. Evolution of /nli, for tiie potential V = Vox'^e'^^ . Comparison with Fig. [3] 
shows that /nl is totally dominated by /nl2- 

5. Discussion 

Non-linearities are now routinely extracted from all-sky observations of the microwave 
background anisotropy. Our purpose in this paper has been to propose a new technique 
with which to predict the observable signal. Present data already give interesting 
constraints on the skewness parameter /nl, and over the next several years we expect 
that the Planck survey satellite will make these constraints very stringent. It is even 
possible that higher-order moments, such as the kurtosis parameter g^i^ [61] will become 
better constrained [65] . To meet the need of the observational community for comparison 
with theory, reliable estimates of these non-linear quantities will be necessary for various 
models of early-universe physics. 

A survey of the literature suggests that the 'conventional' 6N method, originally 
introduced by Lyth & Rodriguez, remains the method of choice for analytical study 
of non-Gaussianity. In comparison, our proposed moment transport method exhibits 
several clear differences. First, the conventional method functions best when we base 
the SN expansion on a flat hypersurface immediately after horizon exit. In our method, 
we make the opposite choice and move the fiat hypersurface as close as possible to 
the time of observation. After this, the role of the 6N formula is to provide no more 
than the non-linear gauge transformation between field fluctuations and the curvature 
perturbation. We substitute the method of moment transport to evolve the distribution 
of field fluctuations between horizon exit and observation. 

Second, in integrating the transport equation one uses an expansion of the velocity 
field such as the one given in Eqs. (!59 l) -( !60l) . This expansion is refreshed at each step 
of integration, so the result is related to conventional perturbative calculations in a 
very similar way to renormalization-group improved perturbation theory f66]. In this 
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interpretation, derivatives of Ui play the role of couplings. At a given order, m, in the 
moment hierarchy, the equations for lower-order moments function as renormalization 
group equations for the couplings at level-m, resumming potentially large terms before 
they spoil perturbation theory. This property is shared with any formalism such as 6N 
which is non-perturbative in time evolution, but may be an advantage in comparison 
with perturbative methods. We also note that although 6N is non-perturbative as a 
point of principle, practical implementations are frequently perturbative. For example, 
the method of Vernizzi & Wands [I9] and Battefeld & Easther [20] depends on the 
existence of quantities which are conserved only to leading order in eN, and can lose 
accuracy after ~ e-foldings. 

Numerical calculations confirm that our method gives results in excellent agreement 
with existing techniques. As a by-product of our analysis, we note that the large 
non-gaussianities which have recently been observed in sum- and product-separable 
potentials [211 122] are dominated by non-linearities from the second-order part of the 
gauge transformation from 6(pi to (. The contribution from intrinsic non-linearities 
of the field fluctuations, measured by the skewnesses aijk, is negligible. In such cases 
one can obtain a useful formula for /nl by approximating the field distribution as an 
exact Gaussian. The non-Gaussianity produced in such cases arises from a distortion of 
comoving hypersurfaces with respect to adjacent spatially flat hypersurfaces. 

Our new method joins many well-established techniques for estimating non- 
Gaussian properties of the curvature perturbation. In our experience, these techniques 
give comparable estimates of /nl, but they do not exactly agree. Each method invokes 
different assumptions, such as the neglect of gradients or the degree to which time 
dependence can be accommodated. The mutual scatter between different methods can 
be attributed to the theory error inherent in any estimate of /nl- The comparison 
presented in §1] shows that while all of these methods slightly disagree, the moment 
transport method gives good agreement with other established methods. 
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